#######################################
#### 06_estimateCJR_model_brazil.R ####
#######################################

rm (list=ls())


print (version)

# platform       aarch64-apple-darwin20      
# arch           aarch64                     
# os             darwin20                    
# system         aarch64, darwin20           
# status                                     
# major          4                           
# minor          1.2                         
# year           2021                        
# month          11                          
# day            01                          
# svn rev        81115                       
# language       R                           
# version.string R version 4.1.2 (2021-11-01)
# nickname       Bird Hippie            



###Vectorization for SIMPLE CJR MODEL  ####

library (runjags)
library (coda)
library (reshape2)
library (random)
library (scales)

#### Call data #### 
savePathRep <- "/all_data/"
# load (paste0(savePathRep, "brazil_CJR_auxiliary_data.Rda"))
load (paste0(savePathRep, "brazil_initial_values.Rda"))
load (paste0(savePathRep, "brazil_rollcall_long_format_mean.Rdata"))

randomseed <- 503006
potential.fix.bill<- c(25, 41, 42,71, 240, 241)



##########################
#### SIMPLE CJR MODEL ####
##########################
CJR.only.v <- "model { 
   for (i in 1:n.obs) {
       y[i] ~ dbern(pi[i])
       probit(pi[i]) <- beta[vote[i]]*theta[dep[i]] + eta[vote[i]]*delta[dep[i]]  - alpha[vote[i]]   
   } 
 # PRIORS
for (j in 1:n.item){ alpha[j] ~ dnorm(0,1) }   
# Beta (discrimination, dimension 1)
for(j in 1:(fix.bill[1]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[1]] ~ dnorm(0,1) 
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[2]]  ~ dnorm(0,1) #
for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[3]] ~ dnorm(0,1)
for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[4]] <- -1.3
for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[5]] <- 1.6
for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { beta[j]  ~ dnorm(0,1) }
beta[fix.bill[6]] <- -3
for(j in (fix.bill[6]+1):n.item) { beta[j]  ~ dnorm(0,1) }
# Eta (discrimination, dimension 2)
for(j in 1:(fix.bill[1]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[1]] <- -2
for(j in (fix.bill[1]+1):(fix.bill[2]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[2]] ~ dnorm(0,1)
for(j in (fix.bill[2]+1):(fix.bill[3]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[3]] ~ dnorm(0,1)
for(j in (fix.bill[3]+1):(fix.bill[4]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[4]]  <- 2 
for(j in (fix.bill[4]+1):(fix.bill[5]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[5]] <- -1.2#
for(j in (fix.bill[5]+1):(fix.bill[6]-1)) { eta[j]  ~ dnorm(0,1) }
eta[fix.bill[6]] ~ dnorm(0,1) #
for(j in (fix.bill[6]+1):n.item) { eta[j]  ~ dnorm(0,1) }
# ideal points
for(i in 1:n.legs)  { theta[i] ~ dnorm(0,1) }
for(i in 1:n.legs)  { delta[i] ~ dnorm(0,1) }
}"

y1 <- molten.roll.call$rc
y2 <- y1
for (i in 1:length(y1)){
  if (is.na(y1[i])==FALSE)
  {if(y1[i]==-5| y1[i]==-7|y1[i]== -9)
    y1[i]<- NA}
  
}

cjr.data.vector <- dump.format(list(y=y1
                                    , fix.bill=potential.fix.bill 
                                    , n.legs=max(molten.party.member$final.leg.no)
                                    , n.item=max(as.numeric(molten.roll.call$vote))
                                    , n.obs=nrow(molten.roll.call)
                                    , vote=as.numeric(molten.roll.call$vote)
                                    , dep=molten.party.member$final.leg.no
))

cjr.parameters = c("theta","delta","alpha","beta","eta","deviance")

set.seed (randomseed)
eta1 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
beta1 = initials$beta + runif(length(initials$beta), -0.2, 0.2)
eta1[potential.fix.bill[4:5]] <- c(NA, NA)
eta1[potential.fix.bill[1]] <- NA
beta1[potential.fix.bill[4]] <- NA
beta1[potential.fix.bill[5:6]] <- c(NA, NA)
CJR.inits.1 <- function() {
  dump.format(
    list(
      theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
      delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
      alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
      eta = eta1,
      beta = beta1
      ,'.RNG.name'="base::Wichmann-Hill"
      ,'.RNG.seed'= randomseed+2
    ))
}
set.seed (randomseed)

eta2 = initials$eta + runif(length(initials$eta), -0.2, 0.2)
beta2 = initials$beta + runif(length(initials$beta), -0.2, 0.2)
eta2[potential.fix.bill[4:5]] <- c(NA, NA)
eta2[potential.fix.bill[1]] <- NA
beta2[potential.fix.bill[4]] <- NA
beta2[potential.fix.bill[5:6]] <- c(NA, NA)
CJR.inits.2 <- function() {
  dump.format(
    list(
      theta = initials$theta + runif(length(initials$theta), -0.2, 0.2),
      delta = initials$delta + runif(length(initials$delta), -0.2, 0.2),
      alpha = initials$alpha + runif(length(initials$alpha), -0.2, 0.2),
      eta = eta2,
      beta = beta2
      ,'.RNG.name'="base::Wichmann-Hill"
      ,'.RNG.seed'= randomseed+3
    ))
}


set.seed(randomseed+1)
CJR.model <- run.jags(
  model=CJR.only.v,  # add .small for small sample
  monitor=cjr.parameters,
  method="parallel",
  n.chains=2,
  data=cjr.data.vector,
  inits=list (CJR.inits.1(), CJR.inits.2()),
  thin=30, burnin=4000, sample=100,
  check.conv=FALSE, plots=FALSE)
date()

chainsCJR <- mcmc.list(list (CJR.model$mcmc[[1]], CJR.model$mcmc[[2]]))  
gelman.CJR <- gelman.diag (chainsCJR, multivariate=F) 
max(gelman.CJR[[1]][,1], na.rm=T)

# Uncomment next line to save file
# save (chainsCJR, file = "brazil_CJR_auxiliary_output.Rda")
